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Abstract — A peak-seeking control method is presented which 
utilizes a linear time-varying Kalman filter. Performance func- 
tion coordinate and magnitude measurements are used by the 
Kalman filter to estimate the gradient and Hessian of the 
performance function. The gradient and Hessian are used to 
command the system toward a local extremum. The method 
is naturally applied to multiple-input multiple-output systems. 
Applications of this technique to a single-input single-output 
example and a two-input one-output example are presented. 

I. Introduction 

Peak-seeking control, also referred to as extremum- 
seeking control, is concerned with on-line optimization 
across an unknown performance function of some physical 
system or process. Given noisy coordinate and magnitude 
measurements of the performance function, it attempts to 
move a plant to a local extremum. It has been studied in re- 
lation to many applications, such as formation flight for drag 
reduction [1], bio-reactors [2], and axial-flow compressor [3]. 
Peak-seeking control gained much attention in the 1950s and 
1960s [4], [5], [6], [7] but interest in this field declined in the 
1970s and 1980s. It has since had a resurgence of attention 
in the past 15 years [2], [8], [9]. 

In general, peak-seeking control methods can be divided 
into three categories: classical-gradient methods, paramet- 
ric methods, and nonlinear methods. The classical-gradient 
methods estimate the performance function gradient using 
classical-control techniques and move the system accord- 
ingly [10]. Parametric methods parameterize the performance 
function to estimate the extremum position and move the 
system accordingly [1]. Nonlinear methods utilize nonlinear 
control techniques such as adaptive control and numerical 
optimization with line searches to estimate the gradient and 
move the system towards the extremum [11], [9]. 

The earliest work used classical-gradient methods partly 
due to ease of implementation in analog devices and the 
lack of inexpensive digital computers. The classical-gradient 
methods have recently been further developed in work by 
Ariyur [12], Krstic [13], and others who have developed a 
general design technique and proved local convergence of 
this approach [10]. The parametric methods are a newer ap- 
proach. Some examples of work in this area are Najson [14], 
Chichka [1], Popovic [8]. Examples of nonlinear approaches 
can be found in [11] and [9]. 

This paper discusses a parametric method of peak-seeking 
control. It uses a linear time- varying Kalman filter to estimate 
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the gradient and Hessian of the performance function and a 
Newton-method to drive the system toward an extremum. 
The paper is organized as follows. Section 2 provides an 
overview of the method. Section 3 discusses the Kalman 
filter design used to estimate the gradient and Hessian of the 
performance function. Section 4 presents a one-dimensional 
example and section 5 presents a two-dimensional example. 
Section 6 concludes the paper. 

II. Overview of Peak-Seeking Control Method 

The goal of peak-seeking control is to move the position, 
x(t ), of a plant, P, to the extremum of a performance func- 
tion /(x(t)). Figure 1 displays this peak-seeking scheme’s 
overall interconnection. Here, x(t ) is one of a number 
of states of P. Performance function magnitude may also 
affect the plant; in the figure, this is indicated by a dashed 
line. At a given iteration k , the difference between the 
current plant position, Xk, and the previous position, Xk-i, is 
calculated. The difference between the current performance 
function magnitude, /(£&), and the previous, f(xk- 1 ), is 
also calculated. These differences are used by a linear time- 
varying Kalman filter (KF) to estimate the current gradient, 
bk, and Hessian, M&. The gradient and Hessian are combined 
to form a position command, x c = bkM^ 1 , to drive the 
plant toward the performance function extremum. A filter 
(Filt) smoothes and scales the command to avoid providing 
P large step commands which can create unwanted actuator 
movements. A persistent excitation signal (PE) is added to 
the smoothed position command to ensure observability of 
the performance function. In practice, an initial position 
command initiates movement of the system. The plant is 
assumed to be stabilized by an inner-loop control. 



The method uses position measurements directly to es- 
timate the gradient and Hessian. This requires the use of 
a linear time- varying Kalman filter. Use of a Kalman filter 
allows designers to utilize an array of Kalman filter design 









techniques. This provides flexibility to the scheme and allows 
it to be applied to many problems. 


III. Kalman Filter Design 


Performance function gradient and Hessian are estimated 
using a linear time- varying Kalman filter whose states consist 
of elements at the current position. This is accomplished 
by parameterizing the performance function using first- and 
second-order terms of a Taylor series expansion. 

Consider the Taylor series expansion of a performance 
function f{pc) about Xk 

f(x) « f(x k ) + bl k (x-x k ) 

+ ^(x - x k ) T M Xk (x - x k ) + o(x - x k ) 

where b Xk is the gradient at Xk, M Xk is the Hessian at Xk, 
and o(-) represents higher order terms. Evaluating (1) at Xk-i 
and rearranging yields 

A/fe = b k Ax k + ^Ax k M Xk Ax k (2) 

where Ax k = x k _ x - x k and Af k = /(a; fc _i) - f{x k ). 

The higher-order terms have been dropped by assuming the 
performance function is adequately modeled as a quadratic 
function at any particular position. 

For simplicity, we restrict ourselves to a two-dimensional 
case. Denote the positions in the two dimensions at time k 
as x\ k and X 2 k . Denote the corresponding gradients as b\ k 
and b 2 k and the corresponding Hessian as 


M k 


Further denote 


^llfc M 12 k 

Ml2 k M 2 2 k _ 


Ax lk = x lk _ 1 - x lk 
Ax 2k = X 2 k _ 1 - X 2k . 


Equation (2) is then written as 
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Equation (3) implies that a parameter estimation technique 
may be used to estimate the gradient and Hessian. Since the 
gradient and Hessian may change with x, and measurements 
of Axi k , Ax 2k and Afk may be noisy, a Kalman filter is 
an appropriate choice of estimator. The Kalman filter states 
are chosen to be ^ 
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This allows the measurement equation of a linear time- 
varying Kalman filter to take the form 


A/fe = -fffeCfe + v k (4) 


where 
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and Vk represent a zero-mean Gaussian white-noise process 
with variance Vfe. 

The gradient and Hessian are modeled as a Brownian noise 
process since they may change in an unknown manner with 
x. The Kalman filter process equation is, therefore, given by 


Ofe+i = Kk + w k 


( 5 ) 


where I is a 5 x 5 identity matrix and Wk represents a zero- 
mean Gaussian white-noise process with variance Wk- 
The linear time-varying Kalman filter is implemented with 
the following equations: 


Pk + 1 — Pk + Wk 

Cfe+1 = Ck + P k H^V k -\Af k - H k Q k ) 
P k = {P k 1 + HlV k 1 H k )- 1 


(6a) 

(6b) 

(6c) 


where P is the state covariance matrix, P the predicted 
state covariance matrix, and ( the state vector. The values of 
Wk and Vk are used as tuning parameters. Typically, initial 
guesses of Wk and Vk are based on previously-obtained 
measurements of the noise or on a noise model. A trial- 
and-error process is then used to adjust the values in order 
to improve the estimates. Detailed derivations of the linear 
time- varying Kalman filter can be found in [15] and [16]. 

The Kalman filter may be expanded to include N mea- 
surements at each iteration k. Define 


A f k ,n = f(x lk ,x 2k ) -f(x i k _ n ,x 2k _ n ) (7a) 

Ax H,n = X H ~ X H-n ( 7b ) 

Ax 2ktn = x 2k - x 2k _ n (7c) 

and Vk, n as the corresponding process noise. The index n 
takes values between 1 and N. The expansion is implemented 
by modifying the measurement equation (4) as 


A/fc,i 


~Vk, l" 

A/fe, 2 

— P-kCk V 

Vk, 2 

_A fk,N_ 


_ v k,N_ 


Here 




\ Ax V 

Dk, 1 

Ax lh,l 

Ax 2 k ,-L 

II 

h Ax V 

\ Ax V 

Dk, 2 

Ax l k ,2 

Ax 2 k ,2 


-\ Ax V 

\ Ax V 

Dk,N 

a x lkN 

Ax 2 k . N 


and Dk, n = A^ 2 fc . The process equation remains as 

it is shown in equation (5). The Kalman filter is implemented 



as shown in equation (6) with 


0 
0 

Vk,N_ 

The number of measurements is used as a tuning para- 
meter. A larger N increases the observability and tolerance 
to noise by providing a over-determined set of equations. It 
also increases the area of the performance function to which 
the gradient and Hessian are fit. For a performance function 
in which the Hessian changes as a function of position, a 
too-large N may slow convergence. 
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IV. One-Dimensional Example 

A one-dimensional example of the method is presented 
in this section. Consider the problem presented in [14]. A 
continuous linear plant model given by 
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seeks the extremum of a performance function. In this 
example the performance function was chosen to be f(x) = 
(cos(x/8.4) + 1.5) (x/6 — 0.4) in place of the quadratic 
function used in [14]. This performance function provides a 
gradient and Hessian that change as a function of position. 
The performance function magnitude measurements were 
corrupted with Gaussian distributed noise with a standard 
deviation of 0.1. There was no noise on the position mea- 
surements. The system was implemented in a 1.0 Hz fixed- 
step discrete simulation. The Kalman filter operated at 0.25 
Hz. The matrix in equation (8) was selected to have 10 
rows. The other elements of equation (8) were selected to be 
of compatible size. The command filter was set to 1 and the 
persistent excitation, PE in Figure 1, was set to 0. 

An initial command was provided to the plant. As the 
system responded to the command, position and performance 
function magnitude measurements were provided to the 
Kalman filter. Estimates of the gradient and Hessian from 
the Kalman filter were combined, x c = to command 

the plant toward the local extremum. 

As the system approached the extremum, A fk became 
small and was buried in noise, leading to poor estimates. 
Typically, the Hessian estimate suffers more than the gradi- 
ent estimates. To partially compensate for this, the system 
switched between a steepest-descent approach and a Newton 
approach. When the smallest singular value of the Kalman 
filter’s error covariance, a(Pj . c ) < 0.005, a Newton approach 
was used. When a(Pk) > 0.005, a steepest-descent approach 


was used. The switching threshold was used as a tuning 
parameter and selected by trial and error. 

The plant position as a function of time is presented in 
Figure 2. The minimum position is depicted by a solid line 
and is reached by the system in approximately 80 seconds, 
after which noisy estimates cause the plant to deviate from 
the minimum. 

The gradient and Hessian estimates are shown in Figure 3. 
The Kalman filter estimates are represented by dashed lines. 
The true gradient and Hessian are represented by solid lines. 
The system required three measurements before beginning 
estimation, thus, the figures show the first non-zero estimate 
at 12 seconds. It is apparent that the estimations began to 
suffer once the system neared the minimum. If the simulation 
were allowed to execute for a longer time period, position 
would have randomly moved about the minimum. This ex- 
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Fig. 2. Plant position versus time. Dashed line: Plant position. Solid line: 
Performance function minimum. 



Fig. 3. Estimated and true gradient and Hessian versus time. Dashed line: 
Kalman filter estimates. Solid line: True gradient and Hessian. 

ample illustrated application of the method to a simple one- 
input one-ouput problem. The gradient and Hessian estimates 
track the true values well and the system quickly reaches the 
extremum. The next example illustrates the natural extension 
of the method to multiple-input problems. 


V. Two-Dimensional Example 

A two-dimensional example of the scheme is presented 
in this section. Consider the problem of formation flight for 
drag reduction. In this problem, a trailing aircraft attempts 
to fly its wingtip in the wingtip vortex of a leading aircraft 
to reduce the induced drag. Substantial fuel savings can be 
achieved by employing such a technique. 

This problem has been studied in multiple papers. A 
formation of two Dornier Do 28-D1 aircraft has been flown 
under automatic control achieving 15% flight power reduc- 
tion [17], and two F/A-18 aircraft have achieved 20% drag 
reduction and 18% fuel flow reduction [18]. More details of 
formation flight including theoretical derivations of the drag 
reduction achievable in formation flight can be found in [1] 
and [17]. 

In this example, a two-ship formation of transport class 
aircraft is modeled and a peak- seeking control system is 
designed to maximize drag reduction. It is assumed that the 
lead aircraft flies in a straight-and-level path. This allows the 
vortex generated by the lead aircraft to be modeled as static 
maps of induced drag coefficient and rolling moment on the 
trailing aircraft as a function of lateral and vertical relative 
position. The induced drag coefficient map is used as the 
performance function. The magnitude of the rolling moment 
map for any given position is used as a disturbance input to 
the trailing aircraft model. The maps were generated using 
a vortex-lattice method with the trailing aircraft wingtip 
positioned inside the leading aircraft wingtip vortex. For 
each position of the map, the aircraft was first trimmed for 
straight-and-level flight and then the induced drag coefficient 
and rolling moment were calculated. It is assumed that the 
vortex changes little with respect to relative longitudinal 
spacing. 

The trailing aircraft is modeled with an 11 -state, 4-input, 
10 Hz discrete state- space model. The modeled states are 
body-axis vertical, lateral, and longitudinal velocities; roll, 
pitch, and yaw angles; roll, pitch, and yaw rates; and inertial- 
axis lateral and vertical relative positions between aircraft. 
The inputs are elevator deflection, aileron deflection, rudder 
deflection, and thrust. The effects due to a changing induced 
drag coefficient are not modeled. 

Normally distributed random noise with a standard devi- 
ation of 0.001 is superimposed on the induced drag coeffi- 
cient performance function magnitude. In addition, normally 
distributed random noise with a standard deviation of 0.012 
meter is superimposed on the position measurements. 

A. Inner-Loop Control Design 

The primary goal of the inner-loop control law is to 
move the trailing aircraft to track relative vertical and lateral 
position commands between the leading and trailing aircraft. 
The secondary goal is to minimize roll angle to ensure the 
trailing aircraft wing remains in the vortex during lateral 
movement. The third goal is to maintain a constant relative 
longitudinal velocity to prevent the trailing aircraft from 
slowly drifting out of formation. 


In order to meet these goals, an inner-loop control system 
was designed which penalizes roll angle and change in 
longitudinal velocity. It follows relative lateral and vertical 
position commands by commanding elevator, rudder, aileron, 
and thrust. An FQR-tracker control design methodology was 
selected for construction of the control system. The aircraft 
model was augmented with integral states of the lateral 
position error, vertical position error, longitudinal velocity 
command, and roll angle. Controller gains were computed 
by minimizing the standard FQR cost function 

oo 

x t Qx + u T Ru d t (9) 

where Q and R are designer selected weightings on the 
states, x , and inputs, u, respectively. The resulting gains were 
used in the interconnection shown in Figure 4. In the figure, h 



Eig. 4. Inner-loop-control-system interconnection. 

represents vertical relative position, y lateral relative position, 
( j) roll angle, and u longitudinal velocity. Vertical, lateral, and 
longitudinal velocities; roll, pitch, and yaw angles; and roll, 
pitch, and yaw rates are contained in £. Elevator deflection is 
represented by Se, aileron deflection by 5a , rudder deflection 
by 5r, and thrust by 5t- Control gains on the aircraft states 
are represented by TQ, and control gains on the errors by 
K e . The subscript c on the loop inputs indicates a command 
to the system. 

B. Kalman Filter Design 

A Kalman filter was designed as discussed in section III 
above. It was chosen to iterate at 0.1 Hz to allow the aircraft 
to travel some distance between iterations. Measurements 
were taken at 10 Hz in between the iterations and averaged 
to form /*. )n , x\ kn , and ^ 2 fc , n of (7). The Kalman filter rate 
was used as a tuning parameter and selected by trial and 
error. The matrix in equation (8) was selected to have 
15 rows (TV = 15). The other elements of equation (8) were 
selected to be of compatible size. 

A persistent excitation was chosen as a 3 rad/sec 0.7 
meter sinusoidal signal superimposed on the relative-position 
command. This allowed N values of fk, n , x\ k n , and ^ 2 fc , n 
to be distributed over a full excitation period. 

The command filter was chosen to be an 10 Hz integrator, 
Filt = 0.1 / (2 — 1). This resulted in a ramping position 
command in place of the step command generated by the 
0.1 Hz Kalman filter estimates. In addition, the system was 
implemented to switch between a steepest-descent and a 







Newton approach. When the Hessian estimate was positive 
definite and the minimum singular value of the Kalman filter 
error covariance a(P/ c ) < 6 a Newton approach was used. 
Gradient descent was used when this was not the case. The 
switching threshold was used as a tuning parameter and 
selected by trial and error. 



C. Simulation Results 

The system was tested by simulation. In the simulation, 
the aircraft was initially positioned to the left and above the 
leading aircraft right wingtip vortex core. Figure 5 depicts the 
path the aircraft followed during the peak-seeking simulation. 
The contours represent induced drag coefficient. The system 
was initially commanded to trace a 0.7 meter radius circle to 
generate initial gradient and Hessian estimates. It was then 
allowed to move toward the minimum. The system primarily 
moved orthogonally through the contours of the plot as it 
moved to the minimal location. The system reached the local 
minimum in 300 seconds. 




Fig. 6. Gradient estimates and true gradient. Solid line: truth. Dots: 
Estimate. 
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Fig. 7. Hessian estimates and true Hessian. Solid line: Truth. Dots: 
Estimate. 


Fig. 5. Aircraft path along induced drag coefficient performance function. 

Figures 6 and 7 show the gradient and Hessian estimates 
as a function of time. The solid lines in the figures represent 
the true gradient and Hessian while the dashed lines represent 
the estimates at each Kalman filter iteration. The gradient 
estimate approximates the true gradient well over the length 
of the entire simulation, as shown in Figure 6; however, the 
Hessian estimate illustrated in Figure 7 is less accurate. 

The error between position commands and the aircraft 
response is depicted in Figure 8. The error never exceeds 
0.8 meters. The oscillatory behavior is due to the excitation 
command. 

The aircraft Euler angles are depicted in Figure 9. The 
aircraft angles stay within reasonable values, never exceeding 
5 degrees. The high-frequency oscillatory appearance of the 
angles is due to non- smooth commands being provided to the 
aircraft. Improvement to the command filter could minimize 
these oscillations. The slower-period oscillations are due to 
the excitation. The roll and yaw angles share the task of 
moving the aircraft laterally. By changing the weightings 
contained in R of the inner-loop control design (9), surface 
movements can be tuned to use more roll or yaw angle. 


The aircraft surface deflections are displayed in Figure 
10. As with the Euler angles, the high-frequency oscillation 
is due to the non- smooth commands to the aircraft and the 
slower oscillation is due to the excitation. Aileron deflection 
goes to 10 degrees and rudder deflection goes to 5 degrees 
when the aircraft is tracing the initially-commanded circle. 
The simulation ends with all surface deflections except 
aileron near 0. Aileron deflection remains at 3 degrees be- 
cause the aileron continues to counteract the vortex-induced 
rolling moment. 

This example has demonstrated the application of the 
method to a two-input one-output problem. The method can 
also be applied to problems with larger numbers of inputs; 
however, the number of estimation parameters grows with 
the number of inputs, m, as YhLi i + m - The estimation 
problem then demands a larger number of measurments, N. 

VI. Conclusion 

A parametric peak- seeking control technique has been 
presented. It utilizes a linear time varying Kalman filter 
to directly use coordinate and magnitude measurements 
of the performance function to estimate the gradient and 
Hessian. The Kalman filter is developed from a Taylor 




0 


200 300 400 500 


y — 



200 300 400 

Time (sec) 


Fig. 8. Vertical and lateral position errors. 


Fig. 10. Aircraft surface deflections versus time. 
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Fig. 9. Aircraft Euler angles versus time. 


series approximation of the performance function about the 
current coordinate measurement. It uses the difference be- 
tween subsequent measurements of position and performance 
function magnitude. The estimates are used in a Newton 
method optimization approach to drive the plant toward the 
performance function extremum. Both a one-dimensional 
and a two-dimensional example are provided to illustrate 
the technique. Both demonstrate good results. In the one- 
dimensional example, the scheme reaches the minimum in 
80 seconds and displays good gradient and Hessian esti- 
mates while using noisy performance function magnitude 
measurements. Application to a two-dimensional problem 
is exhibited with the problem of formation flight for drag 
reduction. In this example, the minimum is reached in 300 
seconds and exhibits good gradient and Hessian estimates. 
This example used noisy performance function magnitude 
and coordinate measurements. The flexibility of the linear 
time- varying Kalman filter allows applicability of the method 
to a large number of problems. 
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Example Application 


•! Formation flight for drag reduction 

-iGoal: Position trailing aircraft to achieve optimal 
drag reduction. 

-lExtremum position unknown 


-!Can measure relative position, magnitude of drag 
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Assumptions 


-iPerformance function has single extremum 
-iMeasureable coordinates and magnitude 
-iGausian distributed noise 
— iPlant is stable and controllable 

•! Inner loop control design treated as separate problem 
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Kalman Filter construction 
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Kalman Filter Measurement E 



Taylor Series expansion of PF 
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Improves observability 
Improves tolerance to noise 
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Kalman Filter State E 


Gradient and Hessian change according to 
Brownian motion 
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Formation Flight Application 


Transport Class Aircraft 
Constraints: 

•! Smooth response by aircraft 

— ! Minimize roll angle (use rudder for lateral positioning 
— ISlow and small persistent excitation 

•! Noisy measurements 





Inner Loop Control 


! LQR-tracker design 


/‘•X 


x Qx H - u F\ H clt 


7 


0 


h 






Kalman Filter Design 


•! Induced Drag Coefficient 
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Kalman Filter Design 


•! Iterates at 0.1 Hz 

— iProvide change in magnitude in between 
measurements 

•! KF using 15 measurements { H k^ Rl5x5 ) 

-iTradeoff between improving noise characteristics and 
slower convergence 

•! PE=3 rad/sec 0.7 m sinusoid 

— iMinimizes passenger discomfort while providing 
observability 

• j Fill - 0.1 j {z - 1) 


Vertical Position (m) 


Peak Seeking Results 
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Gradient Estimation 




Hessian Estimation 



Aircraft response 


•! Roll angle kept under 4 degrees 
•! Yaw angle reaches 5 degrees 
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Aircraft Response 




Aileron deflection reaches 9 degrees 
Rudder deflection reaches 10 degrees 
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